<!DOCTYPE html>

<html xmlns="http://www.w3.org/1999/xhtml">

<head>

<meta charset="utf-8" />
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="pandoc" />
<meta http-equiv="X-UA-Compatible" content="IE=EDGE" />

<meta name="viewport" content="width=device-width, initial-scale=1">

<meta name="author" content="Jared Huling" />

<meta name="date" content="2020-01-14" />

<title>Multi-category Treatments with personalized</title>



<style type="text/css">code{white-space: pre;}</style>
<style type="text/css" data-origin="pandoc">
a.sourceLine { display: inline-block; line-height: 1.25; }
a.sourceLine { pointer-events: none; color: inherit; text-decoration: inherit; }
a.sourceLine:empty { height: 1.2em; }
.sourceCode { overflow: visible; }
code.sourceCode { white-space: pre; position: relative; }
div.sourceCode { margin: 1em 0; }
pre.sourceCode { margin: 0; }
@media screen {
div.sourceCode { overflow: auto; }
}
@media print {
code.sourceCode { white-space: pre-wrap; }
a.sourceLine { text-indent: -1em; padding-left: 1em; }
}
pre.numberSource a.sourceLine
  { position: relative; left: -4em; }
pre.numberSource a.sourceLine::before
  { content: attr(data-line-number);
    position: relative; left: -1em; text-align: right; vertical-align: baseline;
    border: none; pointer-events: all; display: inline-block;
    -webkit-touch-callout: none; -webkit-user-select: none;
    -khtml-user-select: none; -moz-user-select: none;
    -ms-user-select: none; user-select: none;
    padding: 0 4px; width: 4em;
    color: #aaaaaa;
  }
pre.numberSource { margin-left: 3em; border-left: 1px solid #aaaaaa;  padding-left: 4px; }
div.sourceCode
  {  }
@media screen {
a.sourceLine::before { text-decoration: underline; }
}
code span.al { color: #ff0000; font-weight: bold; } /* Alert */
code span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code span.at { color: #7d9029; } /* Attribute */
code span.bn { color: #40a070; } /* BaseN */
code span.bu { } /* BuiltIn */
code span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code span.ch { color: #4070a0; } /* Char */
code span.cn { color: #880000; } /* Constant */
code span.co { color: #60a0b0; font-style: italic; } /* Comment */
code span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code span.do { color: #ba2121; font-style: italic; } /* Documentation */
code span.dt { color: #902000; } /* DataType */
code span.dv { color: #40a070; } /* DecVal */
code span.er { color: #ff0000; font-weight: bold; } /* Error */
code span.ex { } /* Extension */
code span.fl { color: #40a070; } /* Float */
code span.fu { color: #06287e; } /* Function */
code span.im { } /* Import */
code span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
code span.kw { color: #007020; font-weight: bold; } /* Keyword */
code span.op { color: #666666; } /* Operator */
code span.ot { color: #007020; } /* Other */
code span.pp { color: #bc7a00; } /* Preprocessor */
code span.sc { color: #4070a0; } /* SpecialChar */
code span.ss { color: #bb6688; } /* SpecialString */
code span.st { color: #4070a0; } /* String */
code span.va { color: #19177c; } /* Variable */
code span.vs { color: #4070a0; } /* VerbatimString */
code span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */

</style>
<script>
// apply pandoc div.sourceCode style to pre.sourceCode instead
(function() {
  var sheets = document.styleSheets;
  for (var i = 0; i < sheets.length; i++) {
    if (sheets[i].ownerNode.dataset["origin"] !== "pandoc") continue;
    try { var rules = sheets[i].cssRules; } catch (e) { continue; }
    for (var j = 0; j < rules.length; j++) {
      var rule = rules[j];
      // check if there is a div.sourceCode rule
      if (rule.type !== rule.STYLE_RULE || rule.selectorText !== "div.sourceCode") continue;
      var style = rule.style.cssText;
      // check if color or background-color is set
      if (rule.style.color === '' && rule.style.backgroundColor === '') continue;
      // replace div.sourceCode by a pre.sourceCode rule
      sheets[i].deleteRule(j);
      sheets[i].insertRule('pre.sourceCode{' + style + '}', j);
    }
  }
})();
</script>



<style type="text/css">body {
background-color: #fff;
margin: 1em auto;
max-width: 700px;
overflow: visible;
padding-left: 2em;
padding-right: 2em;
font-family: "Open Sans", "Helvetica Neue", Helvetica, Arial, sans-serif;
font-size: 14px;
line-height: 1.35;
}
#header {
text-align: center;
}
#TOC {
clear: both;
margin: 0 0 10px 10px;
padding: 4px;
width: 400px;
border: 1px solid #CCCCCC;
border-radius: 5px;
background-color: #f6f6f6;
font-size: 13px;
line-height: 1.3;
}
#TOC .toctitle {
font-weight: bold;
font-size: 15px;
margin-left: 5px;
}
#TOC ul {
padding-left: 40px;
margin-left: -1.5em;
margin-top: 5px;
margin-bottom: 5px;
}
#TOC ul ul {
margin-left: -2em;
}
#TOC li {
line-height: 16px;
}
table {
margin: 1em auto;
border-width: 1px;
border-color: #DDDDDD;
border-style: outset;
border-collapse: collapse;
}
table th {
border-width: 2px;
padding: 5px;
border-style: inset;
}
table td {
border-width: 1px;
border-style: inset;
line-height: 18px;
padding: 5px 5px;
}
table, table th, table td {
border-left-style: none;
border-right-style: none;
}
table thead, table tr.even {
background-color: #f7f7f7;
}
p {
margin: 0.5em 0;
}
blockquote {
background-color: #f6f6f6;
padding: 0.25em 0.75em;
}
hr {
border-style: solid;
border: none;
border-top: 1px solid #777;
margin: 28px 0;
}
dl {
margin-left: 0;
}
dl dd {
margin-bottom: 13px;
margin-left: 13px;
}
dl dt {
font-weight: bold;
}
ul {
margin-top: 0;
}
ul li {
list-style: circle outside;
}
ul ul {
margin-bottom: 0;
}
pre, code {
background-color: #f7f7f7;
border-radius: 3px;
color: #333;
white-space: pre-wrap; 
}
pre {
border-radius: 3px;
margin: 5px 0px 10px 0px;
padding: 10px;
}
pre:not([class]) {
background-color: #f7f7f7;
}
code {
font-family: Consolas, Monaco, 'Courier New', monospace;
font-size: 85%;
}
p > code, li > code {
padding: 2px 0px;
}
div.figure {
text-align: center;
}
img {
background-color: #FFFFFF;
padding: 2px;
border: 1px solid #DDDDDD;
border-radius: 3px;
border: 1px solid #CCCCCC;
margin: 0 5px;
}
h1 {
margin-top: 0;
font-size: 35px;
line-height: 40px;
}
h2 {
border-bottom: 4px solid #f7f7f7;
padding-top: 10px;
padding-bottom: 2px;
font-size: 145%;
}
h3 {
border-bottom: 2px solid #f7f7f7;
padding-top: 10px;
font-size: 120%;
}
h4 {
border-bottom: 1px solid #f7f7f7;
margin-left: 8px;
font-size: 105%;
}
h5, h6 {
border-bottom: 1px solid #ccc;
font-size: 105%;
}
a {
color: #0033dd;
text-decoration: none;
}
a:hover {
color: #6666ff; }
a:visited {
color: #800080; }
a:visited:hover {
color: #BB00BB; }
a[href^="http:"] {
text-decoration: underline; }
a[href^="https:"] {
text-decoration: underline; }

code > span.kw { color: #555; font-weight: bold; } 
code > span.dt { color: #902000; } 
code > span.dv { color: #40a070; } 
code > span.bn { color: #d14; } 
code > span.fl { color: #d14; } 
code > span.ch { color: #d14; } 
code > span.st { color: #d14; } 
code > span.co { color: #888888; font-style: italic; } 
code > span.ot { color: #007020; } 
code > span.al { color: #ff0000; font-weight: bold; } 
code > span.fu { color: #900; font-weight: bold; }  code > span.er { color: #a61717; background-color: #e3d2d2; } 
</style>




</head>

<body>




<h1 class="title toc-ignore">Multi-category Treatments with <code>personalized</code></h1>
<h4 class="author">Jared Huling</h4>
<h4 class="date">2020-01-14</h4>


<div id="TOC">
<ul>
<li><a href="#example-with-multi-category-treatments"><span class="toc-section-number">1</span> Example with multi-category treatments</a></li>
<li><a href="#more-details-on-propensity-scores-for-multi-category-treatments"><span class="toc-section-number">2</span> More details on propensity scores for multi-category treatments</a></li>
</ul>
</div>

<div id="example-with-multi-category-treatments" class="section level1">
<h1><span class="header-section-number">1</span> Example with multi-category treatments</h1>
<p>To demonstrate how to use the <code>personalized</code> package with multi-category treatments, we simulate data with three treatments. The treatment assignments will be based on covariates and hence mimic an observational setting with no unmeasured confounders.</p>
<div class="sourceCode" id="cb1"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb1-1" data-line-number="1"><span class="kw">library</span>(personalized)</a></code></pre></div>
<div class="sourceCode" id="cb2"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb2-1" data-line-number="1"><span class="kw">set.seed</span>(<span class="dv">123</span>)</a>
<a class="sourceLine" id="cb2-2" data-line-number="2">n.obs  &lt;-<span class="st"> </span><span class="dv">250</span></a>
<a class="sourceLine" id="cb2-3" data-line-number="3">n.vars &lt;-<span class="st"> </span><span class="dv">15</span></a>
<a class="sourceLine" id="cb2-4" data-line-number="4">x &lt;-<span class="st"> </span><span class="kw">matrix</span>(<span class="kw">rnorm</span>(n.obs <span class="op">*</span><span class="st"> </span>n.vars, <span class="dt">sd =</span> <span class="dv">3</span>), n.obs, n.vars)</a>
<a class="sourceLine" id="cb2-5" data-line-number="5"></a>
<a class="sourceLine" id="cb2-6" data-line-number="6"><span class="co"># simulated non-randomized treatment with multiple levels</span></a>
<a class="sourceLine" id="cb2-7" data-line-number="7"><span class="co"># based off of a multinomial logistic model</span></a>
<a class="sourceLine" id="cb2-8" data-line-number="8">xbetat_<span class="dv">1</span> &lt;-<span class="st"> </span><span class="fl">0.1</span> <span class="op">+</span><span class="st"> </span><span class="fl">0.5</span> <span class="op">*</span><span class="st"> </span>x[,<span class="dv">1</span>]  <span class="op">-</span><span class="st"> </span><span class="fl">0.25</span> <span class="op">*</span><span class="st"> </span>x[,<span class="dv">5</span>]</a>
<a class="sourceLine" id="cb2-9" data-line-number="9">xbetat_<span class="dv">2</span> &lt;-<span class="st"> </span><span class="fl">0.1</span> <span class="op">-</span><span class="st"> </span><span class="fl">0.5</span> <span class="op">*</span><span class="st"> </span>x[,<span class="dv">11</span>] <span class="op">+</span><span class="st"> </span><span class="fl">0.25</span> <span class="op">*</span><span class="st"> </span>x[,<span class="dv">15</span>]</a>
<a class="sourceLine" id="cb2-10" data-line-number="10">trt.<span class="fl">1.</span>prob &lt;-<span class="st"> </span><span class="kw">exp</span>(xbetat_<span class="dv">1</span>) <span class="op">/</span><span class="st"> </span>(<span class="dv">1</span> <span class="op">+</span><span class="st"> </span><span class="kw">exp</span>(xbetat_<span class="dv">1</span>) <span class="op">+</span><span class="st"> </span><span class="kw">exp</span>(xbetat_<span class="dv">2</span>))</a>
<a class="sourceLine" id="cb2-11" data-line-number="11">trt.<span class="fl">2.</span>prob &lt;-<span class="st"> </span><span class="kw">exp</span>(xbetat_<span class="dv">2</span>) <span class="op">/</span><span class="st"> </span>(<span class="dv">1</span> <span class="op">+</span><span class="st"> </span><span class="kw">exp</span>(xbetat_<span class="dv">1</span>) <span class="op">+</span><span class="st"> </span><span class="kw">exp</span>(xbetat_<span class="dv">2</span>))</a>
<a class="sourceLine" id="cb2-12" data-line-number="12">trt.<span class="fl">3.</span>prob &lt;-<span class="st"> </span><span class="dv">1</span> <span class="op">-</span><span class="st"> </span>(trt.<span class="fl">1.</span>prob <span class="op">+</span><span class="st"> </span>trt.<span class="fl">2.</span>prob)</a>
<a class="sourceLine" id="cb2-13" data-line-number="13"></a>
<a class="sourceLine" id="cb2-14" data-line-number="14">prob.mat &lt;-<span class="st"> </span><span class="kw">cbind</span>(trt.<span class="fl">1.</span>prob, trt.<span class="fl">2.</span>prob, trt.<span class="fl">3.</span>prob)</a>
<a class="sourceLine" id="cb2-15" data-line-number="15">trt.mat &lt;-<span class="st"> </span><span class="kw">apply</span>(prob.mat, <span class="dv">1</span>, <span class="cf">function</span>(rr) <span class="kw">rmultinom</span>(<span class="dv">1</span>, <span class="dv">1</span>, <span class="dt">prob =</span> rr))</a>
<a class="sourceLine" id="cb2-16" data-line-number="16">trt.num &lt;-<span class="st"> </span><span class="kw">apply</span>(trt.mat, <span class="dv">2</span>, <span class="cf">function</span>(rr) <span class="kw">which</span>(rr <span class="op">==</span><span class="st"> </span><span class="dv">1</span>))</a>
<a class="sourceLine" id="cb2-17" data-line-number="17">trt &lt;-<span class="st"> </span><span class="kw">as.factor</span>(<span class="kw">paste0</span>(<span class="st">&quot;Trt_&quot;</span>, trt.num))</a>
<a class="sourceLine" id="cb2-18" data-line-number="18"></a>
<a class="sourceLine" id="cb2-19" data-line-number="19"><span class="co"># simulate response</span></a>
<a class="sourceLine" id="cb2-20" data-line-number="20"></a>
<a class="sourceLine" id="cb2-21" data-line-number="21"><span class="co"># effect of treatment 1 relative to treatment 3</span></a>
<a class="sourceLine" id="cb2-22" data-line-number="22">delta1 &lt;-<span class="st"> </span><span class="dv">2</span> <span class="op">*</span><span class="st"> </span>(<span class="fl">0.5</span> <span class="op">+</span><span class="st"> </span>x[,<span class="dv">2</span>] <span class="op">-</span><span class="st"> </span><span class="dv">2</span> <span class="op">*</span><span class="st"> </span>x[,<span class="dv">3</span>]  )</a>
<a class="sourceLine" id="cb2-23" data-line-number="23"><span class="co"># effect of treatment 2 relative to treatment 3</span></a>
<a class="sourceLine" id="cb2-24" data-line-number="24">delta2 &lt;-<span class="st"> </span>(<span class="fl">0.5</span> <span class="op">+</span><span class="st"> </span>x[,<span class="dv">6</span>] <span class="op">-</span><span class="st"> </span><span class="dv">2</span> <span class="op">*</span><span class="st"> </span>x[,<span class="dv">5</span>] )</a>
<a class="sourceLine" id="cb2-25" data-line-number="25"></a>
<a class="sourceLine" id="cb2-26" data-line-number="26"><span class="co"># main covariate effects with nonlinearities</span></a>
<a class="sourceLine" id="cb2-27" data-line-number="27">xbeta &lt;-<span class="st"> </span>x[,<span class="dv">1</span>] <span class="op">+</span><span class="st"> </span>x[,<span class="dv">11</span>] <span class="op">-</span><span class="st"> </span><span class="dv">2</span> <span class="op">*</span><span class="st"> </span>x[,<span class="dv">12</span>]<span class="op">^</span><span class="dv">2</span> <span class="op">+</span><span class="st"> </span>x[,<span class="dv">13</span>] <span class="op">+</span><span class="st"> </span></a>
<a class="sourceLine" id="cb2-28" data-line-number="28"><span class="st">    </span><span class="fl">0.5</span> <span class="op">*</span><span class="st"> </span>x[,<span class="dv">15</span>] <span class="op">^</span><span class="st"> </span><span class="dv">2</span> <span class="op">+</span><span class="st"> </span><span class="dv">2</span> <span class="op">*</span><span class="st"> </span>x[,<span class="dv">2</span>] <span class="op">-</span><span class="st"> </span><span class="dv">3</span> <span class="op">*</span><span class="st"> </span>x[,<span class="dv">5</span>]</a>
<a class="sourceLine" id="cb2-29" data-line-number="29"></a>
<a class="sourceLine" id="cb2-30" data-line-number="30"><span class="co"># create entire functional form of E(Y|T,X)</span></a>
<a class="sourceLine" id="cb2-31" data-line-number="31">xbeta &lt;-<span class="st"> </span>xbeta <span class="op">+</span><span class="st"> </span></a>
<a class="sourceLine" id="cb2-32" data-line-number="32"><span class="st">    </span>delta1 <span class="op">*</span><span class="st"> </span>((trt.num <span class="op">==</span><span class="st"> </span><span class="dv">1</span>) <span class="op">-</span><span class="st"> </span>(trt.num <span class="op">==</span><span class="st"> </span><span class="dv">3</span>) ) <span class="op">+</span><span class="st"> </span></a>
<a class="sourceLine" id="cb2-33" data-line-number="33"><span class="st">    </span>delta2 <span class="op">*</span><span class="st"> </span>((trt.num <span class="op">==</span><span class="st"> </span><span class="dv">2</span>) <span class="op">-</span><span class="st"> </span>(trt.num <span class="op">==</span><span class="st"> </span><span class="dv">3</span>) )</a>
<a class="sourceLine" id="cb2-34" data-line-number="34"></a>
<a class="sourceLine" id="cb2-35" data-line-number="35"></a>
<a class="sourceLine" id="cb2-36" data-line-number="36"><span class="co"># simulate continuous outcomes E(Y|T,X)</span></a>
<a class="sourceLine" id="cb2-37" data-line-number="37">y &lt;-<span class="st"> </span>xbeta <span class="op">+</span><span class="st"> </span><span class="kw">rnorm</span>(n.obs, <span class="dt">sd =</span> <span class="dv">2</span>)</a></code></pre></div>
<p>We will use the  version of the treatment status vector in our analysis, however, the integer values vector, i.e. , could be used as well.</p>
<div class="sourceCode" id="cb3"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb3-1" data-line-number="1">trt[<span class="dv">1</span><span class="op">:</span><span class="dv">5</span>]</a></code></pre></div>
<pre><code>## [1] Trt_2 Trt_3 Trt_1 Trt_1 Trt_1
## Levels: Trt_1 Trt_2 Trt_3</code></pre>
<div class="sourceCode" id="cb5"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb5-1" data-line-number="1"><span class="kw">table</span>(trt)</a></code></pre></div>
<pre><code>## trt
## Trt_1 Trt_2 Trt_3 
##    81    92    77</code></pre>
<p>Then we construct a propensity score function that takes covariate information and the treatment statuses as input and generate a matrix of probabilities as output. Each row <span class="math inline">\(i\)</span> of the output matrix represents an observation and each column <span class="math inline">\(j\)</span> is the probability that the <span class="math inline">\(i\)</span>th patient received the <span class="math inline">\(j\)</span>th treatment. The treatment levels are ordered alphabetically (or numerically if the treatment assignment vector is a vector of integers). Our propensity score model in this example will be a multinomial logistic regression model with a lasso penalty for the probability of treatment assignments conditional on covariate information:</p>
<div class="sourceCode" id="cb7"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb7-1" data-line-number="1">propensity.multinom.lasso &lt;-<span class="st"> </span><span class="cf">function</span>(x, trt)</a>
<a class="sourceLine" id="cb7-2" data-line-number="2">{</a>
<a class="sourceLine" id="cb7-3" data-line-number="3">    <span class="cf">if</span> (<span class="op">!</span><span class="kw">is.factor</span>(trt)) trt &lt;-<span class="st"> </span><span class="kw">as.factor</span>(trt)</a>
<a class="sourceLine" id="cb7-4" data-line-number="4">    gfit &lt;-<span class="st"> </span><span class="kw">cv.glmnet</span>(<span class="dt">y =</span> trt, <span class="dt">x =</span> x, <span class="dt">family =</span> <span class="st">&quot;multinomial&quot;</span>)</a>
<a class="sourceLine" id="cb7-5" data-line-number="5"></a>
<a class="sourceLine" id="cb7-6" data-line-number="6">    <span class="co"># predict returns a matrix of probabilities:</span></a>
<a class="sourceLine" id="cb7-7" data-line-number="7">    <span class="co"># one column for each treatment level</span></a>
<a class="sourceLine" id="cb7-8" data-line-number="8">    propens &lt;-<span class="st"> </span><span class="kw">drop</span>(<span class="kw">predict</span>(gfit, <span class="dt">newx =</span> x, </a>
<a class="sourceLine" id="cb7-9" data-line-number="9">        <span class="dt">type =</span> <span class="st">&quot;response&quot;</span>, <span class="dt">s =</span> <span class="st">&quot;lambda.min&quot;</span>))</a>
<a class="sourceLine" id="cb7-10" data-line-number="10"></a>
<a class="sourceLine" id="cb7-11" data-line-number="11">    <span class="co"># return the matrix probability of treatment assignments</span></a>
<a class="sourceLine" id="cb7-12" data-line-number="12">    probs &lt;-<span class="st"> </span>propens[,<span class="kw">match</span>(<span class="kw">levels</span>(trt), <span class="kw">colnames</span>(propens))]</a>
<a class="sourceLine" id="cb7-13" data-line-number="13"></a>
<a class="sourceLine" id="cb7-14" data-line-number="14">    probs</a>
<a class="sourceLine" id="cb7-15" data-line-number="15">}</a></code></pre></div>
<p>An important assumption for the propensity score is that <span class="math inline">\(0 &lt; Pr(T_i = t | \mathbf X) &lt; 1\)</span> for all <span class="math inline">\(\mathbf X\)</span> and <span class="math inline">\(t\)</span>. This assumption, often called the positivity assumption, is impossible to verify. However, in practice validity of the assumption can be assessed via a visualization of the empirical overlap of our estimated propensity scores to determine if there is any evidence of positivity violations. The  function also allows us to visualize the overlap of our propensity scores for multi-category treatment applications. The following code results in the plot shown Figure .</p>
<div class="sourceCode" id="cb8"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb8-1" data-line-number="1"><span class="kw">check.overlap</span>(<span class="dt">x =</span> x, <span class="dt">trt =</span> trt, propensity.multinom.lasso)</a></code></pre></div>
<div class="figure">
<img src="" alt="Propensity score overlap plot for multi-category treatment data." />
<p class="caption">Propensity score overlap plot for multi-category treatment data.</p>
</div>
<p>Each plot above is for a different treatment group, e.g. the plot in the first row of plots is the subset of patients who received treatment 1. There seems to be no obvious evidence against the positivity assumption.</p>
<p>As the outcome is continuous and there is a large number of covariates available for our construction of a benefit score, we will use the squared error loss and a lasso penalty. The model can be fit in the same manner as for the binary treatment setting, however only linear models and the weighting method are available. Here we can also specify the reference treatment (the treatment that the non-reference treatments are compared with by each benefit score).</p>
<div class="sourceCode" id="cb9"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb9-1" data-line-number="1"><span class="kw">set.seed</span>(<span class="dv">123</span>)</a>
<a class="sourceLine" id="cb9-2" data-line-number="2">subgrp.multi &lt;-<span class="st"> </span><span class="kw">fit.subgroup</span>(<span class="dt">x =</span> x, <span class="dt">y =</span> y,</a>
<a class="sourceLine" id="cb9-3" data-line-number="3">    <span class="dt">trt =</span> trt, <span class="dt">propensity.func =</span> propensity.multinom.lasso,</a>
<a class="sourceLine" id="cb9-4" data-line-number="4">    <span class="dt">reference.trt =</span> <span class="st">&quot;Trt_3&quot;</span>,</a>
<a class="sourceLine" id="cb9-5" data-line-number="5">    <span class="dt">loss   =</span> <span class="st">&quot;sq_loss_lasso&quot;</span>)</a>
<a class="sourceLine" id="cb9-6" data-line-number="6"></a>
<a class="sourceLine" id="cb9-7" data-line-number="7"><span class="kw">summary</span>(subgrp.multi)</a></code></pre></div>
<pre><code>## family:    gaussian 
## loss:      sq_loss_lasso 
## method:    weighting 
## cutpoint:  0 
## propensity 
## function:  propensity.func 
## 
## benefit score: f_Trt_1(x): Trt_1 vs Trt_3,  f_Trt_2(x): Trt_2 vs Trt_3 
##                f_Trt_3(x): 0 
## maxval = max(f_Trt_1(x), f_Trt_2(x)) 
## which.max(maxval) = The trt level which maximizes maxval
## Trt recom = which.max(maxval)*I(maxval &gt; c) + Trt_3*I(maxval &lt;= c) where c is 'cutpoint'
## 
## Average Outcomes:
##                Recommended Trt_1 Recommended Trt_2 Recommended Trt_3
## Received Trt_1   4.3602 (n = 22) -16.7813 (n = 53)  -18.0823 (n = 6)
## Received Trt_2 -11.3726 (n = 49)   0.4374 (n = 20)   -19.28 (n = 23)
## Received Trt_3 -21.6418 (n = 41) -23.2637 (n = 28)    4.2641 (n = 8)
## 
## Treatment effects conditional on subgroups:
## Est of E[Y|T=Trt_1,Recom=Trt_1]-E[Y|T=/=Trt_1,Recom=Trt_1] 
##                                          20.1492 (n = 112) 
## Est of E[Y|T=Trt_2,Recom=Trt_2]-E[Y|T=/=Trt_2,Recom=Trt_2] 
##                                          20.8058 (n = 101) 
## Est of E[Y|T=Trt_3,Recom=Trt_3]-E[Y|T=/=Trt_3,Recom=Trt_3] 
##                                           23.1599 (n = 37) 
## 
## NOTE: The above average outcomes are biased estimates of
##       the expected outcomes conditional on subgroups. 
##       Use 'validate.subgroup()' to obtain unbiased estimates.
## 
## ---------------------------------------------------
## 
## Benefit score 1 quantiles (f(X) for Trt_1 vs Trt_3): 
##      0%     25%     50%     75%    100% 
## -53.045  -6.620   3.967  16.523  55.446 
## 
## Benefit score 2 quantiles (f(X) for Trt_2 vs Trt_3): 
##      0%     25%     50%     75%    100% 
## -38.412  -7.542   2.107  12.696  55.281 
## 
## ---------------------------------------------------
## 
## Summary of individual treatment effects: 
## E[Y|T=trt, X] - E[Y|T=Trt_3, X]
## where 'trt' is Trt_1 and Trt_2
## 
##  Trt_1-vs-Trt_3     Trt_2-vs-Trt_3   
##  Min.   :-106.090   Min.   :-76.823  
##  1st Qu.: -13.239   1st Qu.:-15.084  
##  Median :   7.933   Median :  4.213  
##  Mean   :   8.641   Mean   :  5.910  
##  3rd Qu.:  33.045   3rd Qu.: 25.392  
##  Max.   : 110.892   Max.   :110.562  
## 
## ---------------------------------------------------
## 
## 17 out of 30 interactions selected in total by the lasso (cross validation criterion).
## 
## The first estimate is the treatment main effect, which is always selected. 
## Any other variables selected represent treatment-covariate interactions.
## 
## 
## 9 out of 15 variables selected for delta 1 by the lasso (cross validation criterion).
## 
##                                       Trt_1      V1     V2      V3     V4
## Estimates for delta (Trt_1 vs Trt_3) 2.6945 -1.7711 1.8696 -3.6292 0.2637
##                                          V5      V9   V10    V13    V15
## Estimates for delta (Trt_1 vs Trt_3) 2.2749 -2.7176 1.351 0.2342 0.5398
## 
## 8 out of 15 variables selected for delta 2 by the lasso (cross validation criterion).
## 
##                                       Trt_2     V1     V2     V5     V6
## Estimates for delta (Trt_2 vs Trt_3) 2.5083 2.7408 1.3139 -3.055 1.6852
##                                          V8    V11    V12    V13
## Estimates for delta (Trt_2 vs Trt_3) 0.7866 2.0719 1.7789 1.3184</code></pre>
<p>The  function now displays selected variables for each of the two benefit scores and shows the quantiles of each benefit score. We can also plot the empirical observations within the different subgroups using the  function, however now it is slightly more complicated. It appears that the average outcome is higher for those who received the level of the treatment they were recommended than those who received a different treatment than they were recommended. Also note that  returns a  object  and can thus be modified by the user. The below example yields Figure .</p>
<div class="sourceCode" id="cb11"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb11-1" data-line-number="1">pl &lt;-<span class="st"> </span><span class="kw">plot</span>(subgrp.multi)</a>
<a class="sourceLine" id="cb11-2" data-line-number="2">pl <span class="op">+</span><span class="st"> </span><span class="kw">theme</span>(<span class="dt">axis.text.x =</span> <span class="kw">element_text</span>(<span class="dt">angle =</span> <span class="dv">90</span>, <span class="dt">hjust =</span> <span class="dv">1</span>))</a></code></pre></div>
<div class="figure">
<img src="" alt="Individual outcome observations by treatment group and subgroup." />
<p class="caption">Individual outcome observations by treatment group and subgroup.</p>
</div>
<p>To obtain valid estimates of the subgroup-specific treatment effects, we perform the repeated training and testing resample procedure using the  function (here we set the number of replications <code>B=5</code>, but this should be a much larger number in practice) :</p>
<div class="sourceCode" id="cb12"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb12-1" data-line-number="1"><span class="kw">set.seed</span>(<span class="dv">123</span>)</a>
<a class="sourceLine" id="cb12-2" data-line-number="2">validation.multi &lt;-<span class="st"> </span><span class="kw">validate.subgroup</span>(subgrp.multi, </a>
<a class="sourceLine" id="cb12-3" data-line-number="3">    <span class="dt">B =</span> <span class="dv">5</span>,  <span class="co"># specify the number of replications</span></a>
<a class="sourceLine" id="cb12-4" data-line-number="4">    <span class="dt">method =</span> <span class="st">&quot;training_test_replication&quot;</span>,</a>
<a class="sourceLine" id="cb12-5" data-line-number="5">    <span class="dt">train.fraction =</span> <span class="fl">0.5</span>)</a>
<a class="sourceLine" id="cb12-6" data-line-number="6"></a>
<a class="sourceLine" id="cb12-7" data-line-number="7"><span class="kw">print</span>(validation.multi, <span class="dt">digits =</span> <span class="dv">2</span>, <span class="dt">sample.pct =</span> <span class="ot">TRUE</span>)</a></code></pre></div>
<pre><code>## family:  gaussian 
## loss:    sq_loss_lasso 
## method:  weighting 
## 
## validation method:  training_test_replication 
## cutpoint:           0 
## replications:       5 
## 
## benefit score: f_Trt_1(x): Trt_1 vs Trt_3,  f_Trt_2(x): Trt_2 vs Trt_3 
##                f_Trt_3(x): 0 
## maxval = max(f_Trt_1(x), f_Trt_2(x)) 
## which.max(maxval) = The trt level which maximizes maxval
## Trt recom = which.max(maxval)*I(maxval &gt; c) + Trt_3*I(maxval &lt;= c) where c is 'cutpoint'
## 
## Average Test Set Outcomes:
##                        Recommended Trt_1          Recommended Trt_2
## Received Trt_1   2.52 (SE = 8.53, 9.28%)  -9.69 (SE = 6.25, 19.36%)
## Received Trt_2 -6.79 (SE = 4.04, 11.84%)  -8.91 (SE = 9.44, 16.64%)
## Received Trt_3 -30.87 (SE = 8.16, 9.28%) -27.23 (SE = 11.57, 17.6%)
##                         Recommended Trt_3
## Received Trt_1 -41.02 (SE = 17.46, 3.04%)
## Received Trt_2   -17.2 (SE = 8.75, 9.76%)
## Received Trt_3    -7.56 (SE = 2.69, 3.2%)
## 
## Treatment effects conditional on subgroups:
## Est of E[Y|T=Trt_1,Recom=Trt_1]-E[Y|T=/=Trt_1,Recom=Trt_1] 
##                                  21.66 (SE = 13.79, 30.4%) 
## Est of E[Y|T=Trt_2,Recom=Trt_2]-E[Y|T=/=Trt_2,Recom=Trt_2] 
##                                   9.87 (SE = 14.88, 53.6%) 
## Est of E[Y|T=Trt_3,Recom=Trt_3]-E[Y|T=/=Trt_3,Recom=Trt_3] 
##                                     14.72 (SE = 8.38, 16%) 
## 
## Est of 
## E[Y|Trt received = Trt recom] - E[Y|Trt received =/= Trt recom]:                  
## 11.77 (SE = 9.73)</code></pre>
<p>Setting the  argument above to  prints out the average percent of all patients which are in each subgroup (as opposed to the average sample sizes). We can see that about 58% of patients were recommended treatment 1 and among those recommended treatment 1, we expect them to have larger outcomes if they actually receive treatment 1 as opposed to the other treatments. The estimated effects are positive within all three subgroups (meaning those recommended each of the different treatments have a positive benefit from receiving the treatment they are recommended as opposed to receiving any another treatments).</p>
<p>We can visualize the subgroup-specific treatment effects using  as usual with results shown in Figure :</p>
<div class="sourceCode" id="cb14"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb14-1" data-line-number="1">plv &lt;-<span class="st"> </span><span class="kw">plot</span>(validation.multi)</a>
<a class="sourceLine" id="cb14-2" data-line-number="2">plv <span class="op">+</span><span class="st"> </span><span class="kw">theme</span>(<span class="dt">axis.text.x =</span> <span class="kw">element_text</span>(<span class="dt">angle =</span> <span class="dv">90</span>, <span class="dt">hjust =</span> <span class="dv">1</span>))</a></code></pre></div>
<div class="figure">
<img src="" alt="Validation results for multi-category treatment data." />
<p class="caption">Validation results for multi-category treatment data.</p>
</div>
</div>
<div id="more-details-on-propensity-scores-for-multi-category-treatments" class="section level1">
<h1><span class="header-section-number">2</span> More details on propensity scores for multi-category treatments</h1>
<p>For cases with multi-category treatments, the user must specify a propensity function that returns <span class="math inline">\(Pr(T = T_i | \mathbf X= \mathbf x)\)</span> for patient <span class="math inline">\(i\)</span>. In other words, it should return the probability of receiving the treatment that was actually received for each patient. For example:</p>
<div class="sourceCode" id="cb15"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb15-1" data-line-number="1">propensity.func.multinom &lt;-<span class="st"> </span><span class="cf">function</span>(x, trt)</a>
<a class="sourceLine" id="cb15-2" data-line-number="2">{</a>
<a class="sourceLine" id="cb15-3" data-line-number="3">    df &lt;-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">trt =</span> trt, x)</a>
<a class="sourceLine" id="cb15-4" data-line-number="4">    mfit &lt;-<span class="st"> </span>nnet<span class="op">::</span><span class="kw">multinom</span>(trt <span class="op">~</span><span class="st"> </span>. <span class="op">-</span>trt, <span class="dt">data =</span> df)</a>
<a class="sourceLine" id="cb15-5" data-line-number="5">    <span class="co"># predict returns a matrix of probabilities:</span></a>
<a class="sourceLine" id="cb15-6" data-line-number="6">    <span class="co"># one column for each treatment level</span></a>
<a class="sourceLine" id="cb15-7" data-line-number="7">    propens &lt;-<span class="st"> </span>nnet<span class="op">::</span><span class="kw">predict.nnet</span>(mfit, <span class="dt">type =</span> <span class="st">&quot;probs&quot;</span>)</a>
<a class="sourceLine" id="cb15-8" data-line-number="8"></a>
<a class="sourceLine" id="cb15-9" data-line-number="9">    <span class="cf">if</span> (<span class="kw">is.factor</span>(trt))</a>
<a class="sourceLine" id="cb15-10" data-line-number="10">    {</a>
<a class="sourceLine" id="cb15-11" data-line-number="11">        values &lt;-<span class="st"> </span><span class="kw">levels</span>(trt)[trt]</a>
<a class="sourceLine" id="cb15-12" data-line-number="12">    } <span class="cf">else</span> </a>
<a class="sourceLine" id="cb15-13" data-line-number="13">    {</a>
<a class="sourceLine" id="cb15-14" data-line-number="14">        values &lt;-<span class="st"> </span>trt</a>
<a class="sourceLine" id="cb15-15" data-line-number="15">    }</a>
<a class="sourceLine" id="cb15-16" data-line-number="16">    <span class="co"># return the probability corresponding to the</span></a>
<a class="sourceLine" id="cb15-17" data-line-number="17">    <span class="co"># treatment that was observed</span></a>
<a class="sourceLine" id="cb15-18" data-line-number="18">    probs &lt;-<span class="st"> </span>propens[<span class="kw">cbind</span>(<span class="dv">1</span><span class="op">:</span><span class="kw">nrow</span>(propens), </a>
<a class="sourceLine" id="cb15-19" data-line-number="19">        <span class="kw">match</span>(values, <span class="kw">colnames</span>(propens)))]</a>
<a class="sourceLine" id="cb15-20" data-line-number="20">    probs</a>
<a class="sourceLine" id="cb15-21" data-line-number="21">}</a></code></pre></div>
<p>Optionally the user can specify the function to return a matrix of treatment probabilities, however, the columns <em>must</em> be ordered by the levels of . An example of this is the following:</p>
<div class="sourceCode" id="cb16"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb16-1" data-line-number="1">propensity.func.multinom &lt;-<span class="st"> </span><span class="cf">function</span>(x, trt)</a>
<a class="sourceLine" id="cb16-2" data-line-number="2">{</a>
<a class="sourceLine" id="cb16-3" data-line-number="3">    <span class="kw">require</span>(nnet)</a>
<a class="sourceLine" id="cb16-4" data-line-number="4">    df &lt;-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">trt =</span> trt, x)</a>
<a class="sourceLine" id="cb16-5" data-line-number="5">    mfit &lt;-<span class="st"> </span><span class="kw">multinom</span>(trt <span class="op">~</span><span class="st"> </span>. <span class="op">-</span>trt, <span class="dt">data =</span> df)</a>
<a class="sourceLine" id="cb16-6" data-line-number="6">    <span class="co"># predict returns a matrix of probabilities:</span></a>
<a class="sourceLine" id="cb16-7" data-line-number="7">    <span class="co"># one column for each treatment level</span></a>
<a class="sourceLine" id="cb16-8" data-line-number="8">    propens &lt;-<span class="st"> </span><span class="kw">predict</span>(mfit, <span class="dt">type =</span> <span class="st">&quot;probs&quot;</span>)</a>
<a class="sourceLine" id="cb16-9" data-line-number="9"></a>
<a class="sourceLine" id="cb16-10" data-line-number="10">    <span class="cf">if</span> (<span class="kw">is.factor</span>(trt))</a>
<a class="sourceLine" id="cb16-11" data-line-number="11">    {</a>
<a class="sourceLine" id="cb16-12" data-line-number="12">        levels &lt;-<span class="st"> </span><span class="kw">levels</span>(trt)</a>
<a class="sourceLine" id="cb16-13" data-line-number="13">    } <span class="cf">else</span> </a>
<a class="sourceLine" id="cb16-14" data-line-number="14">    {</a>
<a class="sourceLine" id="cb16-15" data-line-number="15">        levels &lt;-<span class="st"> </span><span class="kw">sort</span>(<span class="kw">unique</span>(trt))</a>
<a class="sourceLine" id="cb16-16" data-line-number="16">    }</a>
<a class="sourceLine" id="cb16-17" data-line-number="17">    <span class="co"># return the probability corresponding to the</span></a>
<a class="sourceLine" id="cb16-18" data-line-number="18">    <span class="co"># treatment that was observed</span></a>
<a class="sourceLine" id="cb16-19" data-line-number="19">    probs &lt;-<span class="st"> </span>propens[,<span class="kw">match</span>(levels, <span class="kw">colnames</span>(propens))]</a>
<a class="sourceLine" id="cb16-20" data-line-number="20">    probs</a>
<a class="sourceLine" id="cb16-21" data-line-number="21">}</a></code></pre></div>
</div>



<!-- code folding -->


<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
  (function () {
    var script = document.createElement("script");
    script.type = "text/javascript";
    script.src  = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
    document.getElementsByTagName("head")[0].appendChild(script);
  })();
</script>

</body>
</html>
